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ABSTRACT 



Wide, fragile binary stellar systems are found in the galactic field, and have re- 
cently been noted in the outskirts of expanding star clusters in numerical simulations. 
Energetically soft, with semi-major axes exceeding the initial size of their birth cluster, 
it is puzzling how these binaries are created and preserved. We provide an interpre- 
tation of the formation of these binaries that explains the total number formed and 
their distribution of energies. A population of weakly bound binaries can always be 
found in the cluster, in accordance with statistical detailed balance, limited at the 
soft end only by the current size of the cluster and whatever observational criteria are 
imposed. At any given time, the observed soft binary distribution is predominantly a 
snapshot of a transient population. However, there is a constantly growing population 
of long-lived soft binaries that are removed from the detailed balance cycle due to 
the changing density and velocity dispersion of an expanding cluster. The total num- 
ber of wide binaries that form, and their energy distribution, are insensitive to the 
cluster population; the number is approximately one per cluster. This suggests that a 
population composed of many dissolved small- A clusters will more efficiently populate 
the field with wide binaries than that composed of dissolved large- A clusters. Locally 
such binaries are present at approximately the 2% level; thus the production rate is 
consistent with the field being populated by clusters with a median of a few hundred 
stars rather than a few thousand. 
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INTRODUCTION 

large fraction of stars, across 
multiple systems (e.g. 
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. . , . iDuquennov 

iFischer fc Marcvlll992l ; ISana fc Evanj|2010h . The most re 



mass es, are found 
Mavorl Il99ll ; 



cent survey of our local volume within 25 pc finds al most 
half of nearby solar-type stars (|Raghavan et all l201Ch arc 
multiple. The formation of multiples is presumably tied to 
the star formation process, and reproducing the statistics of 
stellar multiplicity is a necessity for a complete star forma- 
tion theory. However binaries, especially energetically soft 
binaries with binding energies low compared to the mean 
kinetic energy of their environment, are s ubject to modifica- 
tion due to encounters with other stars l|Heggidll975l ; lHillj 
1 19751 ). The picture is thus complicated somewhat by the 
dynamical proces sing that takes place in a cluste red birth 
environment fe.g. lKroupall200ll ; [Parker et alj|2009l l. 
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It is particularly interesting to consider very wide bina- 
ries, with separations comparable to or larger than a typi- 
cal star forming region (a few 10 4 to 10 5 AU). If they are 
dynamically created, it is unclear how such very soft pairs 
could survive a clustered birth; in a long lived cluster the 
soft binaries will be destroyed by repeated encounters with 
their cluster siblings. 

If the binaries form from the fragmentation of a sin- 
gle collapsing core, the upper limit on the binary separation 
should be roughly the Jeans length, typically ~ 10 4 AU, and 
even this is unlikely since it implies the fragmentation of a 
core that is already rotating near its break up velocity before 
collapse. At separations above a few 10 AU, the Galactic 
environment begins to impose an upper l imit on observed 
binar i es via encounters wi th field s tars (|Retterer fc King 
19821 ; I Weinberg et al.|[l987h . GMCs jMallada fc Fernandez 
20011 ). and at larger separations the Galactic tidal field 



I Jiang fc Tremaine] |2010l '). iKouwenhoven et ail I^OIOT ) pro- 
vide a nice summary of current observations and the ap- 
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parent difficulty of producing binaries with very wide sep- 
arations. The bottom line is that there is a window in the 
range of about 10 4 to 10 5 AU where binaries are observed, 
and whose formation is difficult to explain. 

The puzzling production of these wide binaries is the 
focus of this paper. Recently, two numerical studies noted 
the formation of very wide systems in dispersing clusters. 
iMoeckel fc Batel |201Cl ). adva n cing t he hydrodynamic clus- 
ter formation model of iBatd (|20091 ) with an n-body code 
after removing all the gas, found that some systems with 
semi-major axes > 10 4 AU formed in the expanding halo. 
While the hard-soft borderline separation becomes larger as 
the cluster expands, this is not the solution; the wide bina- 
ries that form in simulations and are seen in the field are 
soft even by the more generous standards of their dispersing 
nat al clusters. 

iKouwenhoven et alJ (|201dl ) specifically explored wide 
binary formation in expanding clusters, using n-body models 
of clusters from both virial and super-virial initial condi- 
tions, and from fractal and spherical initial conditions. They 
found in their final binary distribution two peaks in semi- 
major axis; one associated with the hard binaries that form 
at core collapse and drive the expansion of an isolated clus- 
ter, and another at much larger separations that they termed 
the dynamical peak. They proposed that this peak forms as 
stars that are randomly associated in the phase space of a 
Maxwellian distribution find themselves bound as the cluster 
potential becomes less important in the halo of an expanding 
cluster. 

In this paper we offer an alternative explanation for the 
creation of a population of soft yet permanent binaries dur- 
ing the expansion of a cluster. The instantaneous soft binary 
population in a cluster reflects a statistical detailed balance, 
as pairs of stars are perturbed into and out of formally bound 
states. Any single soft binary is typically short-lived, though 
the overall population is approximately static. The expan- 
sion of the cluster and the consequent lowering of the stellar 
density leads to a constantly changing environment in which 
the binaries find themselves. Binaries that are energetically 
soft by local standards at the time of their formation, with 
semi-major axes similar to the interstellar separation, can 
survive if the timescale on which the density drops is suf- 
ficiently short compared to the disruption timescale of the 
binary. These binaries, effectively frozen out of the binary 
creation-destruction cycle, become a population of energet- 
ically soft, yet permanent binaries. 

We begin in section [2] by describing a suite of n- 
body simulations performed to explore this idea, and by 
discussing the properties of the resultant binary population. 
In section [3] we develop our explanation of the freeze-out 
process, augmented by the appendices, and in section U we 
discuss the implications of this work for generating a field 
population of wide binaries before summarising our conclu- 
sions in section O 



2 THE SIMULATIONS DESCRIBED 

2.1 Numerical setup and method 

We simulated isolate d star clusters using the GPU-enabled 
version of NBODY6 |Aarsethll200(il ,r2003), largely following 



some of the setups used in IKouwenhoven et alJ d201Gn. Our 
initia l density structure was a Plummer sphere 1 Plummetl 
limb in virial equilibrium, and we took our masses from a 
Kroupal 120021 ) distribution in the range 0.08-8.0 Mq. We 
ignored any tidal effects and stellar evolution, we have no 
primordial binaries , and we evolved the clu sters using stan- 
dard n-body units |Heggie fc Mathieulll986l ) where the total 
mass, gravitational constant, and virial radius are related by 
M — G = R v — 1. In these units the initial cluster crossing 
time is 2v2j and the cluster energy is -0.25. The absolute 
value of the specific binding energy of a binary, e, is unity at 
the borderline between a hard and a soft binary, with very 
soft binaries having e« 1. 

We modeled clusters with 7V= 256 and 7V= 1024 (lk) 
for 5 x 10 4 time units, and 8196 (8k) stars, for 10 5 time 
units. For clarity of presentation, we focus mainly on the lk 
runs. To gain statistical leverage, we ran 48 runs at 256 and 
lk and 12 runs at 8k; when discussing bulk cluster proper- 
ties we are referring to the median values over all runs, and 
binary population distributions use the total binary popula- 
tion summed over all runs. To identify wide binaries, we first 
find the nearest neighbour of each star. If two stars are mu- 
tually nearest neighbours, we calculate their energy relative 
to their center of mass, and those pairs that have negative 
energy are deemed binaries. 

While we present most of our results in n-body units, 
the scaling of n-body simulations to physical units is 
straightforward. Physical lengths are obtained by simply 
multiplying all length scales by the initial virial radius of 
the cluster. Our choice of a mass function give a mean 
stellar mass m w 0.45 Mq. With mass measured in Solar 
masses and distances in parsecs, characteristic velocities in 
km s are obtained by a w 0.066(iVm/i?„) 1/2 , and one unit 
of time is converted to Myr with the factor 15(A r m/7?^) _1/ ' 2 . 
Choosing either an initial velocity dispersion of a few km 
s _1 or an initial virial radius of the order 0.1 pc are equiva- 
lent, and when we discuss results in physical units we choose 
R v = 0.1 pc as an illustrative case. 



2.2 Cluster expansion 

The clusters undergo relaxation driven global expansion as a 
result of stellar dynamical effects. All the results we describe 
here may not apply to expansion from other processes, such 
as gas removal. The overall structure of the clusters proceeds 
along expected lines, which we describe briefly. In figure[T]we 
show the evolution of the Lagrangian radii for the lk runs. 
Lagrangian radii enclose a fixed mass fraction relative to the 
cluster density center, which is determined via the prescrip- 
tion of lCasertano fc Hutl (|l985t ). We exclude escaping stars 
from the calculation of the Lagrangian radii, although they 
are not removed from the simulation and can contribute the 
binary counts. T he long-term evolution o f these radii is well- 
established (e.g. iBaumgardt et al.ll2002T ). The early evolu- 
tion is marked by contraction of the innermost Lagrangian 
radii, as two-body relaxation conducts kinetic energy out- 
ward from the core, which responds by contracting in a vain 



1 For these purposes, any star outside the 0.9 Lagrangian radius 
with a velocity exceeding twice the escape velocity at its current 
radius is deemed to be an escaper. 
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Figure 1. Lagrangian radii for the lk clusters, shown as the 
median value of the 48 runs; solid lines show the radius enclosing 
the labeled mass fraction. The formation time and location of all 
binaries present at the end of the simulations are shown as dots. 
Soft binaries, with e < 1, are shown in orange, and hard binaries 
with e > 1 are shown in grey. 



attempt to find equilibrium (|Lvnden- Bell fc Eggletonlll980t ). 
The evolutionary timescale is approximately the half-mass 
relaxation time, 



t rh = 0.138- 



(1) 



(Gm) 1 /2i n ( 7 iV)' 

where Rh is the half-mass radius, N the number of stars, 
and m the mean stellar mass. We take as the argument in 
the Coul omb logarithm 7 = 0.0 2, appropriate for our mass 
function (|Giersz fc Heggij|l996l ). Core collapse of a Plum- 
mer sphere with equal masses takes place after about 15 
t r h, although the presence of a mass function accelerates the 
collapse as mass segregation drives the massive star s to the 
core. Once there, the mass-segregation instability (|Spitzerl 
Il969h drives the core to collapse faster than in the equal- 
mass case, and core collap se occurs at less than a relaxation 
time (|Giirkan et al.ll2004r i. 

For our lk systems, the initial relaxation time is t r h,o ~ 
32. Core collapse occurs at about t cc as 27, after which bi- 
naries in the core drive the expansion by heating their sur- 
roundings as they harden. As two-body relaxation conducts 
this heat outwards, the cluster expands in response on the 
relaxation timescale. This expansion is approximately self- 
similar, and the post core-collapse evolution can be reason- 
ably described by 



Tl = Tl,0 1 + X 



trh,0 



(2) 



where Tl is some Lagrangian radius. For the half- mass radius 
itself, one can analytically show that S = 2/3, and while this 
value fits well interior to the half-mass radius, the outer radii 
expand with a slightly steeper power over the length of time 
we simulate. The factor x depends on the cluster mass func- 



tion; [GMes^etjjl] (|2010t ) give % » 0.1(m maa; /m) ' 7 , which 
we use here to find x ~ 1-15. 



2.3 Binary properties 

In figure [I] we plot, in addition to the Lagrangian radii, the 
time of formation and radius in the cluster of all the bi- 
naries that are present at the end of the simulation. While 
we identify the final binaries at the end of the simulation 
(t = 5 X 10 4 for the lk runs), we only plot those binaries 
that formed before t = 2.5 x 10 4 , guaranteeing that they 
are truly long-lived binaries and not a transient population. 
This population of binaries is referred to throughout the pa- 
per as 'permanent'. There are two clear populations; those 
that form in the center of the cluster and that are supplying 
the energy to drive the cluster expansion, and a class of wide 
binaries that form in the outskirts of the cluster. These are 
the soft binaries that primarily concern us here. 

In figure [2] we plot the distribution of semi- major axes 
for the permanent binaries, and the spectrum of binary 
energies for three sets of binaries. The black energy his- 
togram shows the permanent binaries; the medium gray his- 
togram shows all binaries (i.e. bound nearest neighbours) at 
t = 2.5 x 10 4 regardless of their eventual permanence or de- 
struction (the 'instantaneous' binary population); and the 
lightest gray shows the distribution of all bound pairs at 
t = 2.5 x 10 4 regardless of their proximity (the 'statistical' 
binary population). The latter two populations will be dis- 
cussed in the fol lowing section. The final d istributions are 
characterized, as lKouwenhoven et al.l |2010l ) pointed out, by 
a peak of hard binaries at small semi- major axes, and a peak 
at large separations which is the focus of that work and 
this one. The borderline between a hard and soft binary is 
cih-s ~ RvN , and so evolves with time in the same sense 
as the cluster expansion. When plotted as their energy spec- 
trum they are revealed as a distribution that rises smoothly 
towards low binding energy. 

The relation between these three populations can be 
seen in figure [5] The statistical population extends to the 
lowest energies, bounded only by the size of the cluster. The 
instantaneous population is further restricted by the defini- 
tion that the stars need to be nearest neighbours. At late 
times the difference between the instantaneous and perma- 
nent populations is only clear at the softest energies. 



3 EXPLAINING THE SOFT BINARY 
POPULATION 

3.1 Freezing-out of detailed balance 

In a star cluster with a Maxwellian velocity distribution, any 
two stars have a statistical chance of having instantaneously 
negative relative energy. In an infinite, uniform distribution 
of stars, there is a population of soft binaries at all energies, 
maintained by detailed balance as stars are perturbed into 
and out of bound states. While an individual binary is a 
transient object, the statistical distribution is steady; this 
has been known for a long time (|Heggielll975r i. The simple 
statistics that describe the very soft energies become com- 
plicated when one moves closer to the hard-soft border, but 
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Figure 2. Semi-major axis and binding energy distribution of the binary population for the lk runs. The semi-major axis plots shows 
only the permanent binaries. The energy plot shows the statistical population in light gray, the instantaneous population in medium 
gray, and the permanent population in black. These plots are normalized by the total number of runs. For reference, the upper scale on 
the semi- major axis plot gives the values for a cluster with an initial virial radius of 0.1 pc; this can be linearly scaled to any preferred 
cluster size. 



the complete steady-state distribution of binaries in a clus- 
ter with_ajmifonri^taticdensity was thoroughly described 
bv lGoodman fc Hud (|l993h . In that idealized density distri- 
bution, a small fraction of the soft binaries are perturbed 
into hard binaries, and become permanent. The statistically 
constant supply of transient soft binaries provides a source 
of hard binaries, though individual soft binaries tend to be 
short lived. 

In the isolated, modest- N clusters we have simulated 
here, the expansion of the cluster provides another route to- 
ward permanence for a soft binary. As the cluster expands 
and the density drops, a binary that is soft by local stan- 
dards can find itself isolated from the perturbations that 
would destroy it in a static cluster. Provided the density 
drops on a timescale short compared to a binary's destruc- 
tion timescale, it is possible for an individual soft binary to 
drop out of the creation-destruction cycle and become per- 
manent. Given the details of the cluster expansion we can 
calculate this freeze-out binary separation a; TZ , below which 
binaries are likely to survive. We detail this calculation in 
Appendix IA1 

Another necessary criterion for potential permanence is 
that a soft binary have a semi-major axis that is less than 
about the mean interstellar separation a sep at the location of 
its creation. A system with a binary separation larger than 
the local packing of stars is subject to constant perturba- 
tion rather than the discrete perturbations a smaller binary 
experiences. The conditions are ripe for soft, permanent bi- 
nary formation at separations a when a/ rz > a sep ~ a. As 
the cluster's structure evolves, the permanent soft binaries 
that form will be characterized by those regions of the clus- 
ter that meet these conditions. 

To show the relationship between the characteristics of 
the permanent binaries and the evolving cluster, in figure 
[3]we plot the final semi-major axis of the permanent bina- 



formation time / Myr, with Rv = 0.1 pc 




10 1 10 2 10 3 10 4 
formation time / n-body units 

Figure 3. The formation time and semi-major axis of all perma- 
nent binaries from the lk simulations. In gray are energetically 
hard binaries and in orange are the soft binaries, matching figure 
[l] The lines show three quantities at the 0.95 and 0.5 Lagrangian 
radii, with the 0.95 lines above the 0.5 lines. Solid lines show 
the mean interstellar separation; dashed lines show the freeze- 
out radius from equatio rlA8l and dotted lines show the hard-soft 
boundary for stars with the mean stellar mass. 



ries against the time when the pair first became associated. 
We also show a; TZ for the 0.95 and 0.5 Lagrangian radii as 
dashed lines, over the range of time during which the density 
and velocity dispersions are reasonably approximated by a 
power law (figure lATjl . The 0.95 freeze-out radius, which is 
the upper dashed line, is well above any of the binaries that 
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form as well as the interstellar separatiorjf], and is about 3 
orders of magnitude larger than the hard-soft semi-major 
axis at that radius, shown by the upper dotted line. 

We show the freeze-out radius at the 0.5 Lagrangian 
radius, because very few permanent soft binaries form in- 
side the half-mass radius (figure [TJ. This is shown as the 
lower dashed line, only about 1 order of magnitude greater 
than the hard-soft boundary at that radius(the lower dotted 
line). Throughout the evolution of the cluster, a/ rz at the 
half-mass radius is smaller than a 3ep by roughly an order of 
magnitude. The conditions for binaries to drop out of the 
statistical balance are not met, and thus very few perma- 
nent soft binaries form interior to the the half-mass radius. 
In the outer reaches of the cluster a sep <a,f rz for the en- 
tirety of the cluster's evolution, and the semi-major axes of 
the permanent soft binaries formed tracks the evolution of 

3.2 Estimating the number of systems that form 

In this picture, in regions where a scp <af rz most near- 
est neighbours that are instantaneously bound should stay 
bound permanently. In order to estimate the number of soft 
permanent binaries that form, we need to determine the 
fraction of nearest neighbours that are bound. 

We begin by estimating a characteristic interstellar 
separation at the virial radius R v as a 3ep w RyN' 1 ^ 3 . 
The velocity associated with this separation is v se p ~ 
(Gm/dsep) 1 ^ 2 ■ The boundary semi-major axis between a 
hard and soft binary is given by au-s ~ RyN , and is 
related to the velocity dispersion by a w (Gm/ah-s) 1 ^ 2 ■ 

If the nearest neighbour of a star is at something close 
to the interstellar separation, what is the fraction of stars 
in the cluster that have a relative velocity low enough for 
them to be bound? The relative velocities are a Maxwellian 
with dispersion V2a, and we are interested in pairs with 
velocities below y/2v sep . The fraction of bound neighbours 
is the integral over the Maxwellian velocity distribution up 
to that maximum velocity, 

Jbound — ■ _ t — / C X 

2 V 71 " Jo 

where we are integrating over the ratio of the velocity to the 
dispersion. The ratio v sep /a ~ A _1//3 , even for clusters as 
small as A ~ 100 less than about 0.2. In this limit the expo- 
nential term is very close to unity, and the bound fraction 
is well approximated by 




As noted, that ratio of velocities is approximately 
A -1 / 3 , so the bound fraction should be jbound ~ A -1 . A 
point of concern is that we are calculating using a single 
value of these variables for the entire cluster, and taking 
those values at the virial radius rather than the outer re- 
gions of the cluster. However, when the expansion of the 
cluster is driven by relaxation and nearly homologous, the 
local velocity dispersion tracks with the virial velocity, offset 

2 After the power-law behavior of the outer cluster is broken, the 
freeze-out radius becomes even less of a constraint due to the fact 
that the velocity dispersion falls off more slowly with time. 



by a factor of order unity (see figure lA"Tj) . Also note that the 
final result boils down to the ratio of v 3ep to a, which can be 
recast as a ratio of ah-s to a sep . The ratio of these two quan- 
tities is roughly constant over the entirety of the simulations 
at all radii, which can be seen in figure [3] The correction to 
the bound fraction calculation due to variations at different 
radii will be a factor of order unity, not affecting the scaling. 

This argument has implications for the final binary en- 
ergy spectrum and semi-major axis distribution. Provided 
there is some region where a S e P <af rz , as the cluster ex- 
pands and a sep grows larger the permanent binary popula- 
tion is built up from small separations up to large ones. If 
we consider a range of semi-major axes around the critical 
separation at a given time, say a factor of 2, the number 
of bound pairs will be roughly fboundN, or of order unity. 
Thus as the region where soft permanent binaries are cre- 
ated sweeps through larger values, seen in figure [3] as the 
region between the black lines, the number of permanent bi- 
naries formed in any logarithmic separation bin should be of 
order unity. Note that this is independent of the population 
of the cluster. As long as there is a region of the cluster pro- 
ducing soft permanent binaries, details of the cluster density 
structure and population do not matter. 

A flat distribution in the logarithm of the semi- 
major axis means that the energy spectrum should follow 
dAs/de oc e _1 . These scalings are roughly born out by our 
simulations; the best fit to the permanent binary spectrum 
over the soft energies in figure [2] is dAs/de oc e -1,2 , and 
the number of permanent soft binaries in logarithmic bins 
of semi-major axis are within a factor of a few over four 
orders of magnitude. 

The permanent soft population is a fossil record of the 
widest interstellar separation achieved by the regions of the 
cluster from which soft binaries are able to freeze out. The 
final population is independent of the cluster parameters, 
such as the detailed density structure. This is in marked 
contrast to the energy spectra of the statistical and instanta- 
neous populations, which are functions of the finite size and 
geometry of the cluster. The steady-state statistical spec- 
trum derived in a uniform, infinite cluster is d Nb /de oc e _5//2 
at soft energies (|Heggieil 19751 ; iGoodman fc Hutll 19931) . When 
the cluster is finite, this spectrum is truncated at the soft 
end by the size scale of the cluster, which defines the largest 
possible separation of two stars, and thus the lowest energy 
binary attainable. 

When the density field is non-uniform, the slope of the 
statistical energy spectrum changes. For instance, a star sit- 
ting at the center of a cluster with a density profile de- 
creasing outwards will see fewer potential binary partners 
at large separations compared to a uniform distribution (see 
appendix [BJ, which leads to a flattening of the energy spec- 
trum relative to the uniform case. Both of these effects can 
be seen in figure [2] where the statistical population has a 
slope at soft energies more like dA^/de oc e -1,5 , flattening 
further at the softest energies where it is eventually trun- 
cated. The instantaneous spectrum, which is some subset 
of the statistical population, will likewise reflect the current 
density structure of the cluster. The permanent spectrum, 
built up from harder to softer energies over time, should 
generically be dAs/de oc e _1 by the arguments presented in 
this section. 
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10° 10 1 10 2 10 3 10 4 
time / n-body units 

Figure 4. As figure [T] but 

3.3 The 256 and 8k runs 

The overall evolution of the 256 and 8k clusters is similar to 
the lk case. In figure |4] we show the Lagrangian radii and 
formation sites of binaries, analogous to figure [1] Recall that 
the 8k runs are analysed at t — 5 x 10 4 . The expansion of 
the clusters is, as expected, similar to the lk case modulo 
a shift of the time of core collapse related to the N depen- 
dence of the initial relaxation time, and the formation sites 
of permanent soft binaries are similarly concentrated exte- 
rior to the half-mass radius.. Extending the 8k runs several 
orders of magnitude in time would have perhaps been ideal, 
but physical arguments about the relevance of extending the 
runs (see the discussion below) made the computational ex- 
pense of that route seem a questionable investment. 

In figure [S] we show the semi-major axis distributions 
for the 256 and 8k runs. The number of binaries formed 
at large separations, about one per cluster, are consistent 
across all the clusters. The truncation of the distribution for 
the 8k run is because of the shorter length of (logarithmic) 
time those clusters were integrated post core-collapse, and 
because the amount of cluster expansion necessary for the 
interstellar separation to reach a given value increases with 

JVl/3_ 



4 DISCUSSION 

We have been dealing with dimensionless units up to this 
point. When applying these results to a cluster in a galaxy, 
there are two limiting length scales to take into account. 
First, there is the largest binary that we are interested in 
producing. As discussed in the introduction, this is about 
10 AU. Binaries produced above that size scale are not of 
particular interest because their numbers are depleted by 
disruption in the galactic field. Second is the tidal radius of 
the cluster. Binaries are produced in the outer regions of the 
cluster, and as seen in figure Q] these reach radii of the order 




10° 10 1 10 2 10 3 10 4 
time / n-body units 

for the 256 and 8k runs. 

100 n-body units. If this is larger than the tidal radius, these 
binaries will not form, at least not in the relaxation-driven 
expansion scenario we have discussed. 

The first constraint, the upper limit of the binaries 
we are interested in, means that we are only concerned 
with binaries with semi-major axes of a less than a few n- 
body units, thus cutting off the final two bins of figure [2] 
This means that each cluster can be expected to produce 
slightly less than a single binary in the range we are in- 
terested in. The second constraint is that the binary is not 
formed beyond the tidal radius of the cluster. With a tidal 
radius of order 10 pc, our choice of initial virial radius means 
that the cutoff is of order 100 n-body units. Inspecting figure 
[1] we see that binaries continue to be produced closer to the 
half-mass radius through to the end of the simulation inside 
of 100 n-body units, but the production in the outskirts will 
only continue in a truly and unrealistically isolated cluster. 

These two constraints both cut off approximately the 
same binaries from consideration, however. This is to be ex- 
pected, to order of magnitude. Roughly speaking, the clus- 
ter becomes tidally limited as its mean density reaches some 
critical value, and a binary in the field can be thought of as 
limited by a similar density condition. Since the binaries are 
formed at about the interstellar separation, it follows that 
the two constraints affect the same binaries. As is clear from 
figure O after the outer Lagrangian radii are outside a re- 
alistic tidal limit at a few thousand time units, the bulk of 
the binaries with separations less than 10 n-body units have 
already been created, and the story is the same for the 256 
runs. 

Looking to the 8k runs, the situation is somewhat grim- 
mer on the soft-binary production front. The formation site 
of nearly all the binaries are at radii beyond several hundred 
n-body units, flirting with the tidal radius with the same 
initial scaling of R v =0.1 pc. A perhaps more-reasonable 
initial virial radius for this larger cluster of, say, 0.25 pc ex- 
acerbates the problem; more populous clusters need to be in 
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Figure 5. Semi-major axis distribution of the permanent binary populations for the 256 and 8k runs. These plots are normalized by the 
total number of runs. The top scale is in physical units with an initial virial radius of 0.1 pc; these can be linearly scaled to any other 
cluster radius. 



a tidally unrestricted setting in order to expand sufficiently 
to produce the soft binaries we are interested in. It appears 
that when the constraints of the galaxy are included, it is 
the less-populous clusters that are more likely to contribute 
to the population of wide binaries, for two reasons. First, 
with each cluster producing of order a single wide binary, 
10 dispersed clusters of 200 stars will yield more wide field 
binaries than a single ONC of about 2000. Second, the tidal 
limitation restricts the ability of larger clusters to produce 
even their single wide binary. 

In the most recent comprehensive survey of multiplic- 
ity in the Solar neig hbourhood, encompassi ng 454 solar-type 
stars within 25 pc, Ra ghavan et al.l |2010T l find 10 binaries 
with estimated semimajor axes greater than 10 4 AU. Thus 
locally, ~ 2% of field stars are in the very wide systems 
that we have explored here. With approximately one such 
binary being created from each dissolving cluster, a field 
populated by lk clusters would have a 0.1% wide binary 
population, while in a field populated by our 256 star clus- 
ters it would be more like 1%. The statistics of the local 
field are thus broadly consistent with the field being syn- 
thesized by clusters of a few hundred stars, rather than 
a few thousand. T his matches the statistics of clustere d 
star formation sites (|Lada k, Lada|[200^ ; |Porras et alj f2003). 
which appear to hav e a median cluster membership of ~ 300 
|Ad ams et al]|2006l ). We emphasize that the clustered for- 
mation scenario should be thought of in terms of a hierarchi- 
cal structure rather t han as clustered versus isolated modes 
l|Bressert et al.ll2010h . 

This formation process relies upon relaxation-driven ex- 
pansion. Before this point of a cluster's life, there should not 
be any wide binaries. This is cons istent with the lac k of ob- 
served wide binaries in the ONC (|Scallv et alJI 19991 ). and it 
is possible that a small number of wide binaries might be 
formed subsequently during the dissolution of the ONC. It 
has often been noted that the observed binary statistics in 
the ONC are similar to that of the field (thus encouraging 



the notion that the field might be comprised of dissolved 
ONC- type clusters), with the notable exception of the fact 
that the ONC is lacking in wide binaries. Conceivably, there- 
fore, this problem could be solved if the wide binaries are yet 
to form. It is however worth emphasizing that the current 
dynamical status of the ONC is not well known observation- 
ally and that, strictly speaking, our present analysis should 
not be quantitatively applied to clusters that expand from 
super-virial conditions as a result of gas expulsion. 

We have no t included any primordial binaries in these 
simulations. As iKouwenhoven et al.l (|2010T ) point out, the 
effect of primordial binarity depends on the characteristics 
of those binaries. The proper way to set up a primordial 
population is not obvious, especially when considering the 
very early evolution of a cluster where some studies advo- 
cate significant early dynamical processing of the primordial 
population (e.g. lKroupa| [l995; Parker et al. 2009), and oth- 
ers argue that the binary population is largely stable by the 
time gra vitational dynamics be come the dominant physical 
process (jMoeckel &i Batell201Cl ). A full treatment of this is- 
sue is beyond the scope of this paper, but we can draw on 
previous work to speculate on the likely effects. 

IKouwenhoven et al.l (|201dh find that including a primor- 
dial binary population can lead to hierarchical systems con- 
sisting of a primordial binary and a dynamicall-formed very 
wide companion. If the primordial binary is small enough 
that it can approximately be treated as its center of mass, 
i.e. it is energetically hard by local standards, then we expect 
that some fraction of the wide systems produced would be 
members of hierarchical triples (or possibly quadruples). In- 
teractions between more comparably sized binaries will nat- 
urally be more complex; these will be sensitive to the length 
of time the cluster spends pre-core collapse, since this deter- 
mines the amount of processing the primordial binaries may 
undergo. This timescale is highly dependent on the initial 
density structure and mass function of the cluster, and this 
is perhaps an avenue for further work. 
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Regarding the initial cluster structure, we note that 
while we have used a Plummer model as our initial density 
distribution, any standard cluste r setup should yie l d very 
similar post-collapse structures. iBaumgardt et all (|2002| ) 
show by way of example that a f ter co re collapse, a Plum- 
mer sphere and a Wo = 3 King (1966) model develop den- 
sity profiles that are very well matched, provided the time 
of comparison is large compared to the collapse timescale. 
We expect that the outer extent of any expanding cluster 
should satisfy the conditions to produce very wide binaries. 
Due to the nearly self-similar post-collapse evolution, clus- 
ters with extremely high initial densities will evolve through 
the binary-creating stages we simulated, provided the tidal 
limit allows the interstellar separation to reach the requisite 
large value to accommodate the binaries. We further ex- 
pect that subvirial, fractal, or other alternate cluster setups 
should proceed similarly, as the expansion after the densest 
phase of core collapse is driven by the same physics. 

Finally, we note that the number of wide binaries we 
find in these simulations (of order 0.1% for the lk runs) 
is smaller by roughly a n ord er of magnitude than found 
by iKouwenhoven et al.l (|201Ch and iMoeckel fc Bate! (|2010h 
(about 1%). This is due to our insistance on the long-term 
permanence of the binaries, which we check by integrating 
the simulations a factor of two longer than the point at which 
we analyse them. Including all of the instantaneous binaries 
in our analysis boosts the numbers to consistency with those 
previous studies. 



5 CONCLUSIONS 

In an effort to identify the physics behind the formation of 
very wide binaries (semi major axes ~ 10 4 - 10 5 AU), we 
have performed simulations of isolated star clusters under- 
going relaxation-driven expansion after core collapse. Our 
explanation for the formation of these very soft, yet perma- 
nent binaries centers around the effect of the changing clus- 
ter structure on the statistical mechanics that are constantly 
creating and destroying binaries at all energies. Specifically, 
when the local density drops on a timescale short compared 
to the disruption time of a binary, it can be frozen out of 
the statistical equilibrium and become permanent. 

The resultant permanent binary population has some 
interesting characteristics. Namely, the binaries form in the 
outskirts of the cluster, with semi-major axes at approxi- 
mately the local interstellar separation; the final semi-major 
axis distribution is roughly flat in log(a); and the total num- 
ber of binaries formed by any single cluster is of order unity, 
regardless of the population of the cluster. Since approxi- 
mately 2% of local field stars are in binaries with separations 
greater than 10 4 AU, the production mechanism presented 
here could account for that fraction if the field is composed 
of dissolved clusters with a characteristic population of a 
few hundred stars, which is consistent with the statistics of 
clustered star formation sites. 



ACKNOWLEDGMENTS 

Our thanks to Douglas Heggie for his comments on the 
manuscript, to the referee for a helpful report, and to Sverre 
Aarseth for generously sharing his expertise and his CPUs. 



REFERENCES 

Aarseth S. J., 2000, in Gurzadyan V. G., Ruffini R., eds, 
The Chaotic Universe, Proceedings of the Second ICRA 
Network Workshop, Advanced Series in Astrophysics and 
Cosmology, vol.10, Edited by V. G. Gurzadyan and R. 
Ruffini, World Scientific, 2000 NBODY 6: A New Star 
Cluster Simulation Code 

Aarseth S. J., 2003, Gravitational N-Body Simula- 
tions. Gravitational N-Body Simulations, by Sverre 
J. Aarseth, Cambridge, UK: Cambridge University Press, 
November 2003. 

Adams F. C, Proszkow E. M., Fatuzzo M., Myers P. C, 
2006, ApJ, 641, 504 

Bate M. R., 2009, MNRAS, 392, 590 

Baumgardt H., Hut P., Heggie D. C, 2002, MNRAS, 336, 
1069 

Bressert E., Bastian N., Gutermuth R., Megeath S. T., 
Allen L., Evans II N. J., Rebull L. M., Hatchell J., John- 
stone D., Bourke T. L., Cieza L. A., Harvey P. M., Merin 
B., Ray T. P., Tothill N. F. H., 2010, MNRAS, 409, L54 
Casertano S., Hut P., 1985, ApJ, 298, 80 
Duquennoy A., Mayor M., 1991, A&A, 248, 485 
Fischer D. A., Marcy G. W., 1992, ApJ, 396, 178 
Gieles M., Baumgardt H., Heggie D. C, Lamers 
H. J. G. L. M., 2010, MNRAS, 408, L16 
Giersz M., Heggie D. C, 1996, MNRAS, 279, 1037 
Goodman J., Hut P., 1993, ApJ, 403, 271 
Giirkan M. A., Freitag M., Rasio F. A., 2004, ApJ, 604, 
632 

Heggie D. C, 1975, MNRAS, 173, 729 

Heggie D. C, Mathieu R. D., 1986, in P. Hut & 
S. L. W. McMillan ed., The Use of Supercomputers in 
Stellar Dynamics Vol. 267 of Lecture Notes in Physics, 
Berlin Springer Verlag, Standardised Units and Time 
Scales, p. 233 

Hills J. C, 1975, AJ, 80, 809 

Jiang Y., Tremaine S., 2010, MNRAS, 401, 977 

King I. R., 1966, AJ, 71, 64 

Kouwenhoven M. B. N., Goodwin S. P., Parker R. J., 
Davies M. B., Malmberg D., Kroupa P., 2010, MNRAS, 
404, 1835 
Kroupa P., 1995, MNRAS, 277, 1507 
Kroupa P., 2001, MNRAS, 322, 231 
Kroupa P., 2002, Science, 295, 82 
Lada C. J., Lada E. A., 2003, ARA&A, 41, 57 
Lynden-Bell D., Eggleton P. P., 1980, MNRAS, 191, 483 
Mallada E., Fernandez J. A., 2001, in Revista Mexicana 
de Astronomia y Astrofisica Conference Series Vol. 11 of 
Revista Mexicana de Astronomia y Astrofisica, vol. 27, 
Dynamical Evolution of Wide Binaries, p. 27 
Moeckel N., Bate M. R., 2010, MNRAS, 404, 721 
Parker R. J., Goodwin S. P., Kroupa P., Kouwenhoven 

M. B. N., 2009, MNRAS, 397, 1577 
Plummer H. C, 1911, MNRAS, 71, 460 



© 2009 RAS, MNRAS OOO.rflfTol 



Permanent soft binaries in dispersing clusters 9 



Porras A., Christopher M., Allen L., Di Francesco J., 

Megeath S. T., Myers P. C, 2003, AJ, 126, 1916 
Raghavan D., McAlister H. A., Henry T. J., Latham D. W., 

Marcy G. W., Mason B. D., Gies D. R., White R. J., ten 

Brummelaar T. A., 2010, ApJS, 190, 1 
Retterer J. M., King I. R., 1982, ApJ, 254, 214 
Sana H., Evans C. J., 2010. larXiv:1009.4197l 
Scally A., Clarke C, McCaughrean M. J., 1999, MNRAS, 

306, 253 

Spitzer L. J., 1969, ApJ, 158, L139 

Weinberg M. D., Shapiro S. L., Wasserman I., 1987, ApJ, 
312, 367 




time / n-body units time / n-body units 



Figure Al. Lagrangian densities and velocity dispersions for the 
lk clusters, calculated as described in the text. The red lines are 
equations IA1I and IA2I 



APPENDIX A: CALCULATING THE BINARY 
FREEZE-OUT SEPARATION 

Because the cluster is expanding, it is not in the homoge- 
neous s tate from which the stea dy-state distribution is de- 
rived bv lGoodman fc Hud (|l993l ) . If we place ourselves in the 
shoes of a given soft binary, we have some expected lifetime 
that is a function of our ener gy, the local s tellar density, and 
the local velocity dispersion (Heggie 1975). In an expanding 
cluster, this expected lifetime is a constantly changing quan- 
tity, and may in some cases become longer than the cluster 
expansion timescale on which the local density declines. As 
the cluster expands, progressively larger-separation binaries 
are thus frozen-out from the binary creation-destruction cy- 
cle. Binaries that are energetically soft by their local stan- 
dards are then able to survive for long times. 

Given the power-law expansion of the Lagrangian radii 
given by equation [2] we can make predictions of this freeze- 
out process. In what follows, for simplicity we take all stars 
to be the mean mass in the cluster. The density at a La- 
grangian radius (we will call this a 'Lagrangian density') 
that follows equation [2] will evolve as 



for a given semi-major axis and approximated as 



Gma. 



20V2T 
3a 



(A3) 



Collecting constant terms as A, we can write the evolu- 
tion of a number of binaries at that semi-major axis as 



XN a 1 + x- 



t-tr. 
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The solution to this equation for times greater than 
some time t* is 



lnJV (t) = N a (t*) 

t* - tr 



Xtrh,0 
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with the exponents related by a = 38. We empirically find a 
similar relation for the velocity dispersion at a Lagrangian 
radius, so that 



°l = o~l,o 1 + X 



t t C c 
trhfi 



-P 



(A2) 



When calculating the Lagrangian dispersion at a given ra- 
dius, we consider stars inside a shell of width 0.02 around 
that radius. For example, the dispersion at the 0.95 radius 
is calculated using the stars between the 0.94 and 0.96 radii. 
We subtracted off the mean radial velocity of the shell, so 
that the bulk expansion of the cluster does not figure in to 
any calculations of local quantities. The evolution of the La- 
grangian densities and dispersions are shown in figure IA11 
as well as the fits given by equations I Al I and I A2I for the 0.5 
and 0.95 Lagrangian radii. In practice we normalized the fits 
to produce the best match, and we only use them where it is 
visually reasonable (thus the limited range of the fit to the 
0.95 dispersion in the lk run). 

The destru ction rate for a soft binary at a given energy 
is calculated by iHeggid (|l975l ). and can be cast into a rate 



For 1 — a + ft < the solution is, roughly speaking, an 
asymptotic limit at t — >■ oo (the first term in square brackets) 
and a decay to that limit from the initial population (the 
second term). The timescale for decay is something like t r hfi> 
which is short compared to length of the simulation, so the 
behavior of the asymptotic limit appears to be the most 
interesting feature. 

We define the time when a population of binaries at 
some semi-major axis is frozen out as the time when their 
asymptotic limit is e _1 . While this is somewhat arbitrary, 
other reasonable choices (such as one-half surviving) intro- 
duce only a small correction in what follows. Reintroducing 
the numerical factors contained in A, we find that for a La- 
grangian shell with initial density no and velocity dispersion 
ctq the largest semi-major axis that is frozen out evolves as 



frz 



(*) = 



-3 X go(l ~ a + P) 
20v / 2TrnoGmt r h o 



1 + X 



i tcc 
t r h,0 
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This calculation is not a guarantee that a given binary will 
survive, but it provides an estimate of when we might expect 
to see binaries at a given semi-major axis start to become 
permanent at different Lagrangian radii. 
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APPENDIX B: THE STATISTICAL BINARY This paper has been typeset from a TpX/ MftX file prepared 

ENERGY DISTRIBUTION IN A by the author. 

NON-UNIFORM DENSITY FIELD 

At very soft energies, the equilibrium d i stribu tion of bi- 
nary energies found by iGoodman fc HutJ (|l993f) converges 
to a power law with diVs/de oc e~ 5 ^ 2 , a consequence of a 
Maxwellian distribution of velocities and a uniform num- 
ber density l|Heggielll975| y Here we illustrate the effect of a 
non-uniform density distribution, where we assume that the 
typical star sees a number density that is a power law with 
radius, n(r) oc r~ q . 

The absolute value of the specific binding energy of two 
stars of mass m, separation r and relative velocity v is 

2Gm 1 2 1 2 

e = v = x v , 

r 2 2 ' 

where we have introduced x for the gravitational energy of 
the pair. The number density can then be written as n(r) oc 
x q . With a Maxwellian distribution of relative velocities with 
dispersion a, the number of potential binary partners a star 
sees in the range e to e + de is 



N t oc 



dv 

Gm/e 9e 



u 2 exp ( — — - ) r 2 n(r)dr j de, 



with the lower bound on the integral over all radii given 
by the zero- velocity binary pair at energy e. Note that n(r) 
is not the density structure of the cluster, but rather the 
density profile seen by a star, weighted by all the stars in 
the cluster. Our assumption is that this averaged profile is 
fit by a power law over some range. For a Plummer sphere, 
for instance, q ~ f is appropriate for r comparable to the 
length scale of the distribution. 

Converting from radius to x and including the power- 
law density assumption, we have 



N e oc (^j^ i x ~ e ) 1//2 ex P ( 2 2 ) x<1 4 d a; ^ de 
Ve now introduce y = x — e, and rewrite the integ 
N t cc QH y^exp (y + e^^dy) de. 



The exponential term limits the contribution to the in- 
tegral to values of y < a 2 , and since we are concerned with 
very soft energies we have e a 2 . We can consider the be- 
havior of the integrand in two limits. When a 2 ^> y ^> e, the 
integrand is approximately y q ~ 7 ^ 2 . When a 2 3> e ^> y, the 
limiting behavior of the integrand is y 1 ^ 2 . Thus for q <C 7/2, 
the dominant contribution to the integral occurs with val- 
ues y ~ e. At this point the integrand is turning over to the 
y ^> e limit, y q ~ 7 ^ 2 , and the integral is approximately 

iV E oc e 9 ' 5/2 de. 

While this is approximate due to the unpalatable nature 
of the integral, it recovers the uniform density result when 
q — 0. As q increases the peak of the integrand broadens, 
and this treatment becomes troublesome. However, empiri- 
cally a Plummer sphere yields q ~ f, and we feel comfort- 
able estimating the integral in this fashion. With q ~ f , we 
find the approximate diVs/de oc e -1 ' 5 seen in our statistical 
spectrum over a significant range in figure [2] 
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